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Abstract 

Monte-Carlo methods for zero energy quantum scattering are developed. 
Starting from path integral representations for scattering observables, we 
present results of numerical calculations for potential scattering and scatter- 
ing off a schematic ^He nucleus. The convergence properties of Monte-Carlo 
algorithms for scattering systems are analyzed using stochastic differential 
equation as a path sampling method. 
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I. INTRODUCTION 



Quantum scattering problems play an important role in many branches of physics. Al- 
most all our knowledge about composite microscopic systems is obtained from scattering 
experiments. The theoretical analysis of these experiment is particularly difficult if the tar- 
get has internal degrees of freedom which couple strongly to the projectile. Examples of 
this situation are antiproton-nucleus and /^"-nucleus scattering. Since nucleon-antiproton 
or nucleon-i^~ scattering length is of the same order of magnitude as the typical internu- 
cleon distance []T]-§|, one expects that nuclear scattering of these projectiles is dominated by 
multiple scattering correlations involving many nucleons. 

The standard approach to low energy nuclear multiple scattering is the construction of 
effective optical one-body potentials [Q. It has been possible to describe succesfully the 
phenomenology of elastic scattering and inelastic reactions to specific channels [|I],0 as well 
as exotic atoms 0] in the framework of optical potentials. 

One essential shortcoming of the optical model is that it is very restricted from the point 
of view of nuclear dynamics. Phenomenological difficulties in the description of the strong 
annihilation in antiproton-nucleus and K~ nucleus scattering are attributed to missing dy- 
namical properties of simple optical potentials . Recently it has been shown that a model 
which is dynamically richer than the optical potential is able to produce the observed large 
annihilation 0J^. 

This example shows that the development of alternative methods for multiple scattering 
problems is necessary to understand important aspects of nuclear reactions. For many-body 
problems with a discrete spectrum considerable progress has been achieved by using Monte- 
Carlo techniques [pl,|lO|. They have been applied to quantum field theories atomic 
nuclei |jl^,|l3|, quantum liquids [p!^-[T7[| and other problems in condensed matter physics 
18| , p!9|l . However, for scattering problems with a many-body system as target, standard 



Monte-Carlo methods are not applicable. The fact that scattering wave functions are not 
normalizable, and that the spectrum is continuous, requires the development of methods 
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which are specifically designed for scattering problems. 

In this paper, methods are presented, which allow for the calculation of scattering ob- 
servables at vanishing projectile energy. Their convergence behaviour is estimated using 
stochastic differential equations as a path sampling method. A schematic many-body prob- 
lem is discussed in order to demonstrate that the methods developed for potential scattering 
are applicable to more complicated systems. 

This paper is organized as follows: In section |T| three different path integral represen- 
tations of scattering observables at zero energy are derived. In section |T| numerical algo- 
rithms are developed from these path integral representations and applied to scattering off a 
Gaussian potential. The convergence properties for scattering problems are compared with 
systems which posses a bound state. In section |I^ scattering off a schematic "^He nucleus is 
treated. A discussion of the results is given in 0. 

II. SCATTERING OBSERVABLES AND THE PATH INTEGRAL 

A. Path integrals 

In this section we briefly review the path integral formalism following the presentation 



m 



20|. The Monte-Carlo methods developed in this paper are based on a path integral 



representation of the imaginary time propagator: 

p{xf,Xi\P) = {xf\exp{-PH)\xi) 
For a Hamilton operator H which is the sum of kinetic energy and a local potential V{x) 

H = Ho + Vix) = ^ + Vix) 
one finds from the Trotter product formula |^ 

N-l . / ^ X 3/2 . 

p{xf,Xi\p) = hm^ n J d-^^j \27ie) ^^P(~'^[^]) =• / Dxexp{-S[x]) (1) 

7V-1 ^ N-l 

S[x] = So[x] + Sv[x] So[x] = J2 7ri^n+i - Xnf Sv[x] = eV{xn), (2) 

n=0 n=0 



with xo = Xi, xn = Xf and e = jS/N. 

Using the spectral decomposition of the propagator one sees that the ground state wave 
function can be extracted from it in the hmit (3 ^ oo. For a system with a discrete spectrum 
one finds: 



Xi\[3) = exp{-pEn) '^=°° i'o{xf)ijl{xi) exp(-/5£'o) (3) 

n=0 



ipn are energy eigenstates of the Hamiltion operator H and En the corresponding eigenvalues. 
Therefore solving the high dimensional integral (|l]) yields information about the ground state 
of a physical system. Numerical algorithms for calculating the integral (|l]) can be formulated 
using importance sampling techniques [|10| . The convergence rate in f3 in (Q) is controlled by 



the energy gap between the ground state and the first excited state. Excited states decay 
exponentially. This is a typical property of systems with a discrete spectrum. 

The spectrum of a scattering system in absence of a bound state is purely continuous, 
i.e. there is no energy gap. Like in the case of a discrete spectrum the imaginary time 
propagator can be written in terms of a complete set of eigenfunctions of the Hamiltonian 

El 

p{xf,Xi\P) = J d^k{xf\ij^){ipt\xi)exp{-p^) ^=°° (^^^ ^Lo(^i)^fc=o(^/)- (4) 

The states lip'^) are scattering wave functions with projectile energy E = k'^/{2m). They 
are solutions of the Lippmann-Schwinger equation p^ : 



where \k) is a momentum eigenstate. Alternatively we could have used solutions \ip^) with 
ingoing wave boundary conditions 

or an arbitrary linear combination of these functions in equation (H). 
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From equation (^) one sees that for systems with a continuous spectrum the ground 
state can also be found from the large (3 limit of the imaginary time propagator. In nu- 
merical calculations two closely related problems occur: The scattering wave functions are 
not normalizable and due to the lack of an energy gap excited states are only suppressed 
like 1/(3. Therefore long imaginary times are necessary. In a Monte-Carlo algorithm this 
causes large fluctuations which have to be supressed by generating many path samples. In 
order to use Monte-Carlo methods for systems with a continuous spectrum it is particularly 
important to construct estimators for physical observables with a small variance in order 
to avoid long computation times. In the following sections we derive relations between the 
imaginary time propagator and scattering observables which can be used to find estimators 
for these observables. 



B. The cross section at = 



In equation is used together with the asymptotic form of the wave function at 
E = 

^toW™l + ^; a(E = 0)=47^a^ (5) 

to derive a variational principle for the scattering length a and the cross section a. For 
a numerical calculation of the scattering length with Monte-Carlo methods, the following 



relation between scattering length and scattering wave function is more suitable as it 
does not require the limit r ^ oo: 

"^-^^Ito/^l'/to). (6) 

With this expression one can show that the cross section aX E = can be found from a 
matrix element of the imaginary time propagator 

a{E = 0) = lim {2nm[3f''^2p{2nf{k = 0\Ve-'^^V\k = 0). (7) 

This can be expressed as a path integral by using (|lD: 



{2-Kf{k = 0\Ve-^^V\k = 0) = / d^Xfd^XiV{xf)V{xi) H'^ Dxexp{-S[x]) 
This result is generalized to scattering off a target with internal degrees of freedom in section 
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C. The phase shift 

An alternative path integral representation of the scattering length can be found from 
the equation pT| , [23[ 

AZ(,) ^ TVe-- - TVe- ^ g Hl±i ,k.M-P^^, (8) 

which relates the scattering phase shifts 5i{k) of a system with its statistical properties. Like 
(I), this equation holds if there are no bound states of the projectile. Note that both traces 
in the equation (P) are infinite for a system with a purely continuous spectrum. 

Using the Trotter product formula we find again a path integral representation for (§) 

AZ(/3) = j Dxexv{-S\x]) 

S[x] = So[x] - ln(e-^^l"] - 1). (9) 

The paths are periodic in (3, i.e. x^ = Xq. Problems which are related to the logarithm in 
the exponent of will be discussed later. In the limit (3 —>■ oo one finds from the threshold 
behaviour of the scattering phase shifts limk-,oSi{k) = aik'^^'^^ that AZ is proportional to 
the s-wave scattering length a = Oq 

D. The scattering length and the virial theorem 

An alternative method for calculating of the scattering length can be found from the 
scale transformation properties of the transition matrix T. The transition matrix is defined 

as 



T{E)\k) = V\^ljt). (11) 



One can show that 24 



^p{pf\T{Em = -{4~f^\o{x)\4t^) (12) 



P = \Pf \ = \P^\ = (13) 

o{x) = 2V{x) + X ■ VV{x). (14) 

A proof of this formula is given in appendix 0. 

With this equation one can derive a path integral relation for the scattering length a. 
First the spectral decomposition of the imaginary time propagator in terms of the scattering 
wave functions I'lp'^) and lip^) is used: 



J d^x{xf\exp{—^H)\x)o{x){x\exp{—^H)\xi) p^^oo /l67rm\ 



3/2 



(^,-=o|o(x)|V'to)- (15) 



{xf\eM-m\Si) V /5 I 

Because the wave function \ip'l^) and \ipk) identical at /c = the contributions from the 
points Xf and Xi cancel. 

With the help of the Trotter product formula the lefthand side of the equation can be 
rewritten path integral: 



SS:ffDxeM-S[x])o{x{^)) 



Ihs. = ^ ^ =: (o(x(-)). (16) 

Dxexp(-5N) ' ^ ^2^^ 



With equations (^) and (|TID one finally finds for the scattering length: 



327r(2m)V2 ^ 



III. MONTE-CARLO EVALUATION OF THE PATH INTEGRAL 

A. Reference potentials 

In path integral Monte-Carlo methods a path sampling algorithm is used to generate 
paths which are distributed according to a (normalized) probability weight 
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P(x) = "^P^~^[^]^ (18) 
^ ' JDxexp{-S[x])- ^ ^ 



The average of an estimator E{x) over the sampled paths approaches the expectation value 

m 



1 r 
lim — y,E[xi] = / DxP[x]E[x] =: (E) 



in the limit of an infinite number of samples A^^. In the path integral representations in (|^ 
and (|T0|), the scattering observable is not related to an expectation value of a functional in 
path space but to the path integral itself. Therefore we still have to construct estimators 
for scattering observables. Consider the functional 

N-l 

Ey[x] = expe ^ {V{xn) - V{xn)) = exp(S'y[x]) exp(-S'y [x]), (19) 

n=0 

where V is an arbitrary function. The expectation value of Ey is 

^ J Dxexp{-S[x])Ey[x] ^ J Dx exp{-So[x] - Sy[x]) 
^ J Dxexp{-S[x]) J Dxexp{-So[x]- Sv[x])' ^ ' 

The potential V in the numerator is replaced by the potential V . Therefore from (^ one 
sees that in the limit of large imaginary time (3 the expectation value of Ey is the ratio of 
the ground state wave functions of the two potentials p5[] : 

^k=i)y{Xi)^l=Qy{Xf) 

If a solvable reference potential V is chosen which is similar to V the variance of E will be 
small and the Monte- Carlo algorithm converges rapidly. 

Equation (|2l|) can be used to construct an estimator for the ratio of the cross sections 
of V and V . If the potential V has no zeros we can sample paths which are distributed 
according to the exponential of a modified action 

S'[x\ = S[x] - In V{xf) - In V{xi). (22) 

Unlike the path integral (||) the points Xf and Xi are no longer fixed but also considered as 
integration variables. The expectation value of the estimator 
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(^v^) = ' IZ ^f'^'^J^ (21) 



V{xf)V{x.f''^'^^ ^^^^ 
with paths sampled according to S' is the ratio of the cross sections of V and V 

(E'y) = ^. (24) 
av 

If the potential has zeros one can divide the integral over the endpoints Xf and Xi into parts 
where the potential does not chance sign and perform a simulation for each of these parts. 
In this case the ratio of cross sections is the sum of these contributions. 

An analogous estimator can be constructed for the path integral (P). With the same 
arguments as before one sees that the expectation value of 

= "'pi-^^M'-; (26) 

approaches the ratio of the scattering lengths in the limit of large f5 

(E) ^. (26) 



B. Langevin simulation 

Recently Langevin simulation has been studied as a formal approach to quantum me- 
chanics and as a method for practical calculations (see |2^,27] for a review). The basic idea 



of a Langevin simulation is to solve numerically the stochastic differential equation 

^/nit) = -^S[x{t)] + r/„(t), (27) 
where rjnit) is a Gaussian random variable 

{'nn,i) = {r]n,i{t)r]m,j{s)) = 25nm5ij5{t - s) . 

The time variable t is not the physical time which is represented by the index n here but a 
simulation time. The Langevin equation (^) describes the dynamics of a whole path which 
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moves under the influence of a drift force VnS and random fluctuations rfnit). One can show 
that the stochastic process is ergodic: 

J DxP[x]E[x] = hm ^ dtE[x{t)] (28) 

with P deflned as in (0) and E arbitrary. 

The convergence properties of the Langevin equation are related to the Fokker-Planck 
operator 

which is a differential operator in path space. P{x) is the stationary solution of the Fokker- 
Plank equation 

^$(a;,t)=F$(x,t), (30) 

i.e. FP[x] = 0. The Fokker-Plank equation describes the time evolution of probability 
distributions of paths. If F has a negative discrete spectrum, then $(a;, t) in (^) converges to 
P[x] in the limit t ^ oo if the initial distribution t = 0) is not orthogonal on P. In terms 



of the Langevin equation (p7|) this means that the distribution function of the stochastic 
process Xn{t) in the limit t ^ oo is P[x], independent of the initial distribution. The 
convergence rate, i.e. the typical times T needed in (pS]), is given by the flrst nonvanishing 
eigenvalue of F. 

For a numerical application one discretizes the Langevin time and simulates the stochastic 
process 

dS\x(t)] 

Xn{t + At) = Xn{t) j^At + AuJn{t) (31) 

2At5„^5ij, iit = s 
{AuJn,i{t)AuJmAt)) = \ (32) 

0, else 

We start from an arbitrary initial path af„(t = 0) and iterate ( |3TD with a sequence of random 
numbers Awn{t). The method ( PT| ) corresponds to a simple Euler scheme for the solution of 



ordinary differential equations. Higher order algorithms have been studied recently pq-pO 
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C. An example: scattering off a Gaussian potential 



After having introduced the formal techniques we will now apply these methods to scat- 
tering off a Gaussian potential 

Vix) = Voexpi-^x'). (33) 

As the only relevant parameter of zero energy scattering is 2mVo6^ one can set 6 = 1 and 
m = 1 without loss of generality. 



For the method described in section |1I B|, the action S' in equation (^) is inserted into 



the stochastic differential equation (|27D . As an initial path of we use 

with Xf and Xj randomly chosen inside the potential. Because of the additional part 
\nV{xi) + \iaV{xf) in the action S' (^), the endpoints of the path are subject to a linear 
drift force in the Langevin simulation which moves them back inside the potential region. 
Therefore they will not leave the potential under the influence of the stochastic force Au. 
For numerical integration of the path integral defined in [II Q the action @ is inserted 



into the Langevin equation (pTj). We obtain the stochastic differential equation 

|-W-^-z4l^-^"« .,,.1-.-. (34, 
The logarithm in the action disappears, as only the derivative of the action enters into the 
Langevin equations. For the potential (|33|), which does not change sign, D[x] is always 
positive (y repulsive) or negative (V attractive) as expS'y[x] is either greater or less that 
one. For a general potential D[x], can have zeros. In this case the stochastic process defined 
by (0) is not ergodic. The path space is divided into regions of positive and negative signs 
of D[x]. At the boundaries of this regions the drift force becomes infinite and prevents the 
stochastic process from moving into another region. From the numerical point of view this 
means that several simulations with random initial conditions have to be performed to make 
sure that the whole path space is covered. A more detailed interpretation of the physical 
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meaning of the term in the stochastic differential equation will be given when we 

discuss the convergence properties of the algorithms. For the equation (Q) we use Xn = 
as initial path. 

An estimate of the typical times /3 required in the numerical calculation can be found 
from the first order Born approximation for the matrix element (0): 

. . (4.,„/.3.v(.))^ (l - jlf^) ^ i2.m^.rnV,^ (l - ) (35) 

For the the choice of parameters 6 = 1, m = 1 the typical physical times will be /5 ~ 100 

Since analytical estimates of the step size e are in general hard to obtain, we tested 
discretization errors numerically. For the present work, we used e = 0.5, which resulted in 
errors less than 5 percent. Typically 25000 warmup steps were taken to equilibrate, followed 
by 250000 steps of measurement at At = 0.04. This corresponds to a Langevin time of 
T = 10000. The choice of the Langevin parameters is discussed in section |1IID| where the 
convergence properties of the algorithm are studied in more detail. 
As a reference, we used a square well potential 

\ r> R 

V{r) = . (36) 

l^Vo r<R 

From equation ( |35D we see that for not too strong potentials, the cross section and the 
leading term in 1/(3 is given by the zeroth and second moment of the potential. If the 
parameters R and Vq are chosen such that these moments are identical for both potentials, 
we expect that the ratio of cross sections is close to 1 and depends weakly on (3. This leads 
to an initial guess of the parameters of V 



R = V5b Vo = 3^—Vo. (37) 

We first discuss the /3 dependence of the result for the potential strengths Vq = 0.1, 
Vo = —0.3 and Vq = —0.5. Figure |1] shows the ratio of the cross sections as a function of /3. 
The lines are the analytic result which are obtained by numerically solving the Schroedinger 
equation for both potentials. The data point correspond to the stochastic result. 
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The ratio of the cross section is predicted well for Vq = 0.1 and Vq = —0.3 and does 
not depend on (3 from (3 ^ 10 on. The deviation of the simulation from the exact result 
is typically below 4 percent. For Vq = —0.5 a stronger (3 dependence is found and the 
simulation predicts a result which is approximately 10 percent too large. 

Obviously the reference potential ( P^ which is based on first order Born approximation 
(^) is not a good model for the scattering potential if Vq is large. The Gaussian potential has 
its first bound state at Vq ~ —0.67, whereas the reference potential with the parametrization 
(^) has its first bound state at Vq ~ —0.74. Therefore the reference potential (^) is too 
weak if Vq is close to the bound state value. If we choose the parameters 

R = 2 Vo = ^- (38) 
0.67 32 ^ ^ 

the reference potential also has a bound state at Vq = —0.67. In figure ^ we plot the result 
for this improved reference potential at Vq = —0.5, with all other parameters unchanged. 
The result is now almost (3 independent and deviates by not more than 4 percent from 
the exact calculation (Note the different scales of fig. |I| and fig. |^). This deviation is the 
discretization error in e. 

From the result in figure |^ one sees clearly the importance of choosing a good reference 
potential. With the appropriate choice of V the estimator Ey does not strongly depend 
on j3 and long simulation times can be avoided. In practice one guesses an initial reference 
potential and improves it gradually until the variance of E is minimal and no /3-dependence 
is observed. 

Qualitatively the calculation using the method (|T0|) shows the same (3 dependence as in 
figure and 0. In figure |^ the cross section is plotted for the simulation based on (|^). One 
sees that accurate results can be obtained over a wide range of potential strengths and cross 
sections with the same simulation parameters. 

The central question of every Monte-Carlo calculation is whether the equilibrium dis- 
tribution is sampled correctly and whether enough samples have been chosen to obtain a 
statistically significant result. The answer to this question can be found by studying the 
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autocorrelation time of the algorithm. For Fokker-Planck operators with a discrete spectrum 
one can show that for long times t the correlations of the degrees of freedom decay like 

{Xn{t)Xn{0)) * oT exp(--^), 

lac 

where the autocorrelation time tac = 1/1-^ I is the inverse of the smallest eigenvalue A of the 
Fokker-Planck operator. In the next section we will use this concept to give an estimate of 
the convergence properties of the Langevin simulations. Numerically this question can be 
answered by studying the dependence of the averaged estimator 

e(T) = 1 r dtE'^ixit)] ^ -1 £ E'^[x{U)] (39) 

I JO I\s l^Q 

on the simulation time T. In figure (^ we plot e(T) for Vq = —0.3 for both methods (^ 
and (|l3). One sees strong fluctuations at small T which average out as T increases. From 
T = 4000 on the result is essentially constant. 

From this figure one would estimate an autocorrelation time of the order of magnitude 
of about 1000, i.e. the algorithm produces one independent sample every 25000 paths. This 
has to be compared with typical values for bound state problems in many body physics of 
about 10 — 100 paths for one independent sample (see e.g. 1T1|,[TB| and the references given 
there). These large autocorrelation times are not a problem of the Langevin algorithm, but 
a physical property of scattering systems, as will be seen soon. 

We now discuss the results for method (|1^. For this calculation the initial and final 
point of the path are fixed to Xi = Xf = 0. The stochastic differential equation ( pT]) is 
used for path sampling. In contrast to the previous methods a reference potential is not 
needed as equation ([T7|) directly relates the scattering length a to the expectation value of 
the function o(r). We again plot the result as a function of the simulation time in figure 
(^). The potential strength is Vq = 0.1. Although the potential is rather weak one needs 
T = 100000, i.e. 2.5 x 10^ samples to come close to the analytical result of a = —0.21. 

The slow convergence of (|T7p can be understood from the structure of the estimator o{x) 
and the paths which are generated by the stochastic differential equation. In figure | we plot 
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o(r) and the probability distribution of finding the midpoint of the path at a distance r from 
the origin. Obviously the midpoint of the paths lies outside the potential most of the times, 
i.e. at values of r where o(r) is zero. Only few samples contribute to the expectation value 
(o) in (p!?!). For this reason it is not advisable to calculate the scattering length directly with 
eq. (pri). It can be used together with method (|^ to determine the sign of the scattering 
length. The estimator in equation (|19D yields accurate results for the square of the 

scattering length and the same paths which are used with E'y [x] can be used to measure the 
sign of a via ([17|) . 

D. Estimate of convergence properties 

The numerical results presented in the previous section showed that a rather large number 
of path samples is required in the simulation. Further it was found that the paths are outside 
the range of the potential most of the time. Both observations are intimately related to the 
physical properties of scattering systems. They are not an artifact of the algorithm we use. 

The reason for both problems is that scattering wave functions are not normalizable. 
Consider an attractive potential which vanishes at infinity and is strong enough to have a 
bound ground state at energy Eq. At /3 — ^ oo the typical spatial extension of the paths is 
r ^ 1/ ^/2mEo, as the path integral is dominated by the ground state wave function, which 
decays like ex]i{—y/2mEQr). If the potential is too weak to have a bound state the typical 
extension of the paths diverges as /5 oo, as the zero energy scattering wave function 
(§) does not vanish at r ^ oo. Therefore one expects that the convergence properties of 
the Monte-Carlo simulation becomes worse in the limit of large (3, which means that the 
autocorrelation time strongly depends on (3. 

This can be interpreted intuitively from the properties of the stochastic differential equa- 
tion (^) with the action S' eq. p^) . For a potential V, which is attractive but too weak to 
have a bound state, the potential drift Vn^v is not strong enough to move the paths back 
into the range of the potential. If the endpoints of the paths were not fixed inside the range 
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of the potential by the term \nV{xi) + \nV{xf) in 5", they would move to infinity under 
the infiuence of the stochastic force i]{t). The other points of the paths are kept back in the 
vicinity of the potential by the kinetic energy term 5*0 of the action which causes a linear 
force between neighbouring path points in the stochastic differential equation. For systems 
with at least one bound state the potential energy is strong enough to keep the paths inside 
the range of the potential, even if the endpoints are not fixed. 

One can give a quantitative result for the typical behaviour of the path points by studying 
only the kinetic energy part of the action plus the part of the potential which keeps back 
the endpoints. In this case the stochastic differential equation has the form 



at"'" 



N 



J2 MnmXmit) + f]n{t) 
m=0 



(40) 



771 

M=-(21-A), A:= 

e 



1 - 
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1 

1 1 



(41) 



This equation describes a set of coupled harmonic oscillators. The eigenvalues of the Fokker- 
Planck operator are integer multiples of the eigenvalues of —M [^. Therefore the lowest 
eigenvalue A of M is related to the autocorrelation time tac = 1/A. 

The eigenvalues of M can be found by numerical diagonalization. In figure 0.1 we show 
the autocorrelation time at constant e = 0.5 as a function of /? for different widths h of the 
potential which binds the endpoints. tac diverges like . For our simulation at /? = 100 
we find an autocorrelation time of tac ~ 2000 which agrees with the estimate from plotting 
e(T) in figure ^. 

We also see from figure |^ that the autocorrelation time depends strongly on the width 
h. In the limit 6 — oo the endpoints of the paths are no longer fixed in the range of the 
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potential. In this case the matrix M has one vanishing eigenvalue with the eigenvector 
^min oc (1,1,..., 1). This zero mode corresponds to a displacement of the path. The 
stochastic differential equation is translational invariant in this case. This causes an infinite 
autocorrelation time as the paths spread over a space with an infinite volume. 

Figure (|^,2) shows the same calculation but with a small harmonic oscillator potential 

■1 iV-l 

-e 2^ mu a;„ 

^ n=0 

{uj = 0.1) added to the action. Independent of the conditions on the endpoints we reach a 
constant autocorrelation time at large values of (3, in agreement with our discussion at the 
beginning of this section. 

Therefore two conditions are necessary for the convergence of the Langevin simulation: 
One is a boundary condition on the endpoints of the paths and the other is a finite time 
/3. This guarantuees a discrete spectrum of the Fokker-Planck operator. The situation is 
analogous to quantum mechanics on a finite interval. The boundary conditions on the wave 
function causes the spectrum to be discrete, even for systems which have only scattering 
states in the continuum limit L oo. The energy gap of a system with a bound state 
approaches a nonzero value in this limit, for a system with a purely continuous spectrum it 
vanishes. 

A remarkable point in this result is that a connection exists between the spectrum of 
the Fokker-Planck operator and the Hamiltonian. Although we did not formulate a rigorous 
proof, it is clear from physical arguments that for actions of the form (0) the Fokker-Planck 
operator has a continuous sectrum in the limit /? — > oo, if the Hamilton operator has a 
purely continuous spectrum. 

Let us briefly discuss the method presented in (|T0|) from the point of view of boundary 
conditions on the paths. We sum over periodic paths which are not fixed in space. If a path 
drifts outside the range of the potential, i.e. if Sy in ( p^D approaches zero, the factor 1/D 
increases the potential drift and moves the path back into the range of the potential. This 
is true for repulsive and attractive potentials. Therefore the Fokker-Planck operator of the 
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stochastic differential equation (|3J) has a discrete spectrum for finite /3 and the simulation 
converges. 

Numerically we have seen that the convergence properties of the simulations based on (|^ 
and ( |T0| ) are very similar. The autocorrelation time in both cases is strongly (3 dependent. 
The reason for this can be understood by comparing the potential drifts of (0) with the 
potential drift of e.g. a particle in a harmonic oscillator potential. For the harmonic oscillator 
every point of the path experiences a linear force moving it back to the origin, whereas in 
(^) we have an overall factor which only increases the drift if most of the points lie outside 
the potential. Therefore we still have a diverging autocorrelation time at large (3. 

One can also give an estimate of the maximum Langevin discretization time At from 
the eigenvalues of M. This quantity is determined by the maximum eigenvalue of M as 
convergence of the discretized Langevin equation requires XAt ^ 1 for all eigenvalues A. 



From the Perron- Frobenius theorem for positive matrices ||33| it is found that 4m/e is an 
upper limit for the eigenvalues of M. The step size At therefore has to fulfill the condition 

At < (42) 

For the parameters used here the upper bound is At <^ 0.125. Numerically strong dis- 
cretization errors are found close to this value and instabilities of the simulation above the 
critical value. 



IV. MANY-BODY TARGETS 

In the previous two sections it has been shown that for potential scattering it is possible 
to calculate scattering observables at vanishing energy with Monte-Carlo methods. We now 
give an application to a schematic many-body system in order to show that convergence 
properties are not fundamentally different for targets with internal degrees of freedom. 

Consider a target which consists of a set of A harmonic oscillators with a Hamiltonian 

A 
a=l 
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2M 2 ^'^ , 



(43) 



The projectile interacts with the target via a sum of local two-body potentials 

A 

V{x, gi, . . . , g^) = 51 ^(^ ~ ^^)- 

a=l 

Again the relevant object is the long time limit of the imaginary time propagator. In the 
absence of a projectile target bound state we find analogous to (|) 



(fVI exp(-/3i7)|xg) exp(-/5i?o) f ^ j 7^0,0 (5", g') ^0,0 (x, g), (44) 

where il)ofl{x, q) is a scattering wave function with the projectile at = and the target in 
its ground state |0) as incoming wave. Eq is the ground state energy of the target. We use 
the notation q= (gi, . . . , g^) for the degrees of freedom of the target. 
The elastic cross section of the projectile at vanishing energy is: 

a{E = 0) = lim (27rm^)^/2^e^^»(27r)=^(0, k = 0\Ve-^^V\0, k = 0) (45) 

The matix element can again be written as a path integral 

(27r)^(0, k = 0\Ve-^^V\0, k = 0) = 

A 

J2 J d^^qi£%£xid^XfV{xi - qa,i)%iq.)V{xf - g;',/)$o(g^) x 

DxDqexp{-S[x,q]) . (46) 



a' ,a=l 



^o(g) = (g|0) is the target ground state wave function. 

As the propagator of the harmonic oscillator is known in closed form one can use an 
improved Trotter-product formula and insert for one time step 

(xV|exp(-6if)|fg) ^ Go^°(g',g|e) (^)'^'exp (^-(^(f' - f)^ + 6 ^(f - g^)) j (47) 

with the imaginary time propagator 

- i^^f -p {-j^^if ^ i) -'-^ - Hi)) (48) 

for the harmonic oscillator. Using this path integral means that we have no discretization 
error in the propagation of the target alone, as we summed up the full dynamics of the target 
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already. Unfortunately there are, to our knowledge, no closed expressions for propagators 
of scattering systems in three dimensions. If an exact expression for the dynamics of the 
projectile could be found, only the interplay of the projectile and target degrees of freedom 
would determine the size of the time step e. As the autocorrelation time critically depends 
on the number of time steps for scattering systems, improved short time propagators PT|J5^ 
could be particularly helpful for practical applications. 

The treatment of the endpoints of the paths can be performed in the same way as for 
potential scattering by simulating the modified action 

S'[x, q] = S[x, q] - In $o(£,) - ln$o(£P - In^ ^(^^ " <la,i) ' ^^T.V{x^ - qaj) (49) 

a a 

For a many-body system the integral over the endpoints of the paths is high dimensional 
and can only be handled by stochastic integration. 

In equation (^9]) the logarithm of the ground state wave function is added to the action. 
For the problem discussed here this is a harmonic oscillator wave function which is analyti- 
cally known. At the end of this section we will briefly comment on how this formula can be 
generalized to applications where the ground state wave function of the target is not known. 

As an estimator a functional analogous to (|T^ can be used H 



N-l I A \ 

E\x, g] = exp e ^ ^ V(f„ - q^^a) - Vejjixn) ■ (50) 

n=0 \a=l / 

Here the coupling of projectile and target degrees of freedom is replaced by an effective 
potential V^jj which only contains the projectile degrees of freedom. Unlike the traditional 
approach to nuclear multiple scattering where the many body problem is approximated 
by an effective optical potential, no approximations are made here. The Monte-Carlo method 
yields an exact result and compares it with the effective potential. 

With the same arguments as for potential scattering one can show that measuring the 
estimator (^) with paths sampled with the modified action (|49|) yields in the limit of large 
imaginary time (3 



IJ^oo \/v - w| i^giyfc^o/yyfc^oi vg\n. - w/ ^^^^ 



{k = 0| 






k = 0) 


(0,A; = C 


\V\^o,o){^o,o\V\0,k = 0)' 
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with the folded potential 

Vg{x) = J d'q\%{q)\'J2V{x-q,). (52) 

a=l 

The quantity in the denominator is the cross section of the full many body problem whereas 
the numerator can be found as the solution of a potential scattering problem. If one uses 
the corrected estimator 

El..,]^El.,,]^^^^ (53) 



one again finds that the expectation value of E' is the ratio of the cross sections 

{E') 



a 

There is one major difference as compared to the potential scattering calculation. One can 
always find an estimator Ey with a vanishing variance for potential scattering by setting 
V = V. In other words: we know that there are model potentials which are arbitrarily close 
to the potential which we want to study. This is not true for the estimator (0), because 
correlations of target and projectile degrees of freedom can not be modeled with a one-body 
potential. The variance of E in the simulation will show whether this is possible or not. 

As we are dealing in the following with a system of distinguishable identical particles 
we can further simplify the treatment of the endpoint by decomposing the sum (|45|) into a 
indirect contribution cxj, where the projectile path begins and ends at different nucleons, and 
an indirect contribution ad, where the projectile path begins and ends at the same nucleon 

a{E = 0) = Aaa + A{A - l)a, 

ad = lim (27rm/3)^/2/5e^^"(0, k = 0\V{x- qi)e-'^^V{x - gi)|0, A; = 0) 
ai = lim (27rm/?)^/=^/3e^^o(0, k = (}\V{x - qi)e-'^^V{x - g2)|0, A; = 0). 

The advantage of this procedure as compared to (|^) is that one has already isolated the 
factors A and A{A — 1) correctly. If one simply inserts (^) into the stochastic differential 
equation and starts the calculation with a initial path which begins and ends at the same 
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nucleoli it might take a long time T until the simulation reached a path which begins and 
ends at different nucleons. After splitting the path integral into the direct and indirect 
contribution one can perform two independent simulations with actions 

Sd[x,q] = S[x,q] -ln$o(g) - ln$o(g') - \nV{x-qi) - \nV{x' - g^) 

Si[x,q] = S[x,q] - ln$o(g) - ln$o(g') - In ^(f - gl) - lnV"(f' - ^2)- 

We finally get for the cross section 

(^V,ff }d }i 

One could also specify two independent reference potentials VJ^ and V^jj to define a direct 
and indirect observable. 

In the following this formalism is applied to scattering of a system of four independent 
particles in a harmonic oscillator potential as a simple model for a '^He nucleus. We focus here 
on the convergence properties of the algorithm in the regime of strong multiple scattering 
effect. As interaction between target and projectile the potential ( |33D with b = 0.6 fm is 
used. The mass of projectile and target are set to m = M = 1 GeV ~ 5 fm^^. The 
oscillator frequency is = 0.075 fm~^ which corresponds to a mean square radius of the 
nucleus of 4 fm^ and an energy gap of AE = Ei — Eq = 15 MeV. 

For our calculation we use the folded potential ( ^2|) as an effective potential for the direct 
and indirect observable. If only the ground state of the target contributes as intermediate 
state the folded potential gives the correct cross section. Therefore the deviation of the 
ratio of the cross sections (0) from one is related to the contribution of excited states of 
the target in the Born series. 

The convergence properties of the simulation can again be studied by plotting the average 
of E' ( |53D as an function of the simulation time T analogous to eq.(^). In figure |^ the 
predicted ratio of cross sections for one and four nucleons at /5 = 80 MeV~^ is shown. The 
same discretization was chosen as for potential scattering. 

As expected the result is close to one for the weakly attractive potential Vq = — 10 MeV 
and does not fluctuate strongly. The ground state of the target dominates the Born series 

22 



and the folded potential is a good model. The projectile- nucleus scattering length predicted 
by the Monte-Carlo calculation is ~ 0.75 fm whereas the scattering length off one bound 
nucleon is ai ~ 0.15 fm. This shows that essentially only single scattering occures as 
04 ^ 4ai. 

For the strongly attractive potential Vq = —40 MeV the situation is different. The 
projectile-nucleus scattering length is predicted to be 04 ~ 240 fm and the scattering length 
of a single bound nucleon ai ~ 1 fm. Therefore for this choice of the potential multiple 
scattering effects play a dominant role and the folded potential is no longer a good model 
as excited states of the target become important as intermediate states. 

The statistical fluctuation in this case are larger than for the weak potential Vq = 
—10 MeV. Figure ^ shows, however, that even for Vq = —40 MeV the simulation has a 
stable equilibrium value at T = 10000 fm^. This demonstrates that the method is well 
suited for calculations where multiple scattering effects dominate. 

No qualitative differences are observed in the fluctuations for one and four nucleons. 
This indicates that also simulations with much more than four nucleons are possible, as the 
computational effort scales like the number of target particles. 

Finally we will briefly comment on how the treatment of the endpoints has to be modifled 
if the ground state wave function of the target is not known in a closed form. In this case one 
can use that the Hamiltonian of the target Hf projects a trial state |0) on the true ground 
state 

e-/3i^'|0) ^^=°°e-^i^«|0). 
We can therefore modify equation (|4^) by using the expression 

a{E = 0) = lim (27r)3(27rm/3)^/2^e('^+2/3i)So^0|e-^i^*(/t = 0\Ve-''''V\k = 0)6"'^^^' |0). (55) 

The time f3i has to be chosen large on the scale of the target dynamics. The path integral 
which follows from ( p^D contains projectile paths of length (3 and target paths of length 
f3 + 2f3i. In the simulation we add the logarithm of the trial wave function to the action of 
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the target degrees of freedom in order to fix the endpoints. The folded potential can also be 
found numerically from 

... ft^oo (0|exp(-/3ig,)\/(x,g)exp(-Aif,)|6) 
(0|exp(-2/?i/7,)|0) 

Therefore more complicated Hamiltonians for the target than (^3|) cause no fundamental 
complications. 

V. CONCLUSION 

The basic goal of this work was to show that Monte-Carlo methods can be applied 
succesfully to scattering problems in cases where the target has internal degrees of freedom. 
One possible applications is nuclear multiple scattering, but the methods developed here 
should also be useful in atomic and molecular physics. 

We showed for potential scattering and scattering of a schematic ^He nucleus that despite 
the continuous spectrum of the systems one can perform Monte-Carlo calculations, if one uses 
the concept of a reference potential to reduce statistical fluctuations. For many-body systems 
a one-body reference potential is only a first step in constructing optimal observables. It 
would be useful to find solvable schematic models which take into account the coupling of 
projectile and target degrees of freedom. 

For a phenomenological application of our results a number of questions have to be 
answered: important candidates for calculations in nuclear scattering are antiprotons and 
K~ mesons where the scattering lengths are comparable or larger than typical internuclear 
distances HTj-Q. In both systems annihilation of the projectile plays an important role. This 
is modeled frequently by complex potentials [|^, which are not suitable for Monte-Carlo 



calculations. In order to study these problems one has to find models for annihilation which 
can be implemented in a Monte-Carlo calculation. Another open question is how to deal with 
the spin of nucleons. It is well known that antisymmetrization is very important for nucleons 
as projectiles PH|. A related problem is the implementation of spin-orbit interactions. A 
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problem of all Monte-Carlo algorithms is, that they only converge to the ground state of 
the physical system. Scattering observables can therefore only be calculated in absence of a 
projectile bound state. 

The main advantage of the methods we presented is that one can find numerical solutions 
for a class of many-body problems without approximations on the target dynamics. Ques- 
tions like in-medium effect in nuclear scattering and the reduction to effective potentials can 
be investigated by purely microscopic calculations. 

A further improvement of the methods seems to be possible in several ways: introducing 
improved short time propagators (see e.g. [0,^) will increase the convergence rate of the 



simulations. So far there are very few results on how improved propagator can be constructed 
for scattering systems. One of the great achievements of multiple scattering theory is that 
it is formally possible to separate elementary processes on one nucleon from the multiple 
scattering dynamics [0,^]. A similar formal structure has not yet been developed in a path 
integral framework. 

Another question is whether it is possible to design path sampling algorithms which 
are optimized for scattering problems. As scattering paths move outside the range of the 
potential most of the time it might be possible to find algorithms which use the fact that the 
functional form of the solution in the outside space is known exactly. Langevin simulation 
seems to be a promising method of developing new algorithms because the convergence 
properties of the simulation can be studied in terms of the spectrum of the Fokker-Planck 



operator. For numerical applications, however, hybrid algorithms ||38| , |39| which are based on 
stochastic and deterministic dynamics in path space may be more useful. The applicability 
of these methods to scattering problems is currently studied. 
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APPENDIX A: PROOF OF EQUATION (12) 

Equation can be derived from the scale transformation properties of the Lippmann- 



Schwinger equation ||2J]. The operator 

f/ = g(«P+f)lnA 

generates scale tranformations 

UxU^ = \x UpU^ = \p. 

A 

One can define a rescaled T-matrix and momentum eigenstates: 

{Pi\T{E)\pf) = {\p,\T,iE)\\pf) (A2) 

Tx{E) = UT{E)U^ \Xp) = U\p) 
As the matrix element in equation ( [A^ ) does not depend on A one finds: 

j^{mTx{E)\\pf)=Q (A3) 
From this equation one gets after setting A = 1 



with the definition T' = -^T\ 



^ ^. For this equation we used that T is the solution of the 



operator equation 



Furthermore: 



T{E) = V + V I T{E). 

E — Hq + le 
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In order to prove (|T2D we have to write the derivative with respect to A of the momentum 
eigenstates as derivatives with respect to p. It is easy to show that 



p-^{p,\T{E)\Pf) - m\TxiE)\pf), (A5) 



A=l 

where the operator d/dp only acts on the arguments of the momentum eigenstates. From 
the equations ([A3D , (|A4D , (|A5|) and the definition of the wave functions 

equation ([T2|) follows immediately. 

APPENDIX B: TECHNICAL REMARKS 

All Monte-Carlo simulation in this work were done on an HP9000/720 workstation. A 
simulation at /5 = 100 and e = 0.5 for potential scattering with 250000 path samples requires 
approximately 40 CPU-minutes for method (|^. The code is written in standard FORTRAN 
77. Method (0) needs about the same amount of time. As the convergence properties 
of the many body problem which we studied are comparable to potential scattering, the 
computation time scaled with the number of degrees of freedom. For the ^He calculation 
5 hours of CPU time were needed for one simulation run. Tests with improved random 
number generators and higher order Langevin algorithms showed that an increase of speed 
by a factor of 5 should be possible but we did not study this systematically. 
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FIG. 1. Ratio of the cross sections of the square well potential an the Gaussian potential as a 
function of the physical time /?. The line is the solution of the Schroedinger equation, data points 
are the results of the Langevin simulation. Vq = 0.1 (O), Vq = —0.3 (+) and Vq = —0.5 (□) with 



the reference potential (36 
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FIG. 2. The same calculation as in figure |T| using a modified reference potential ( |38| 
Vo = —0.5. (O): Monte-Carlo result, Line: Schroedinger equation. Note the different scale! 
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cross section of the Gaussian potential 




potential strength VO 
FIG. 3. The cross section of the Gaussian potential as a function of Vq. The hne is the result 

of the Schroedinger equation, data points are the Monte-Carlo simulation (0). 



33 



e (T) for V0=-0 . 3 



0.98 i 



Monte-Carlo (1) 
Monte-Carlo (2) 




2000 



4000 6000 
Langevin time T 



8000 



10000 



FIG. 4. The prediction of the ratio of the cross sections of the Gaussian potential and the 
reference potential as a function of the simulation time T at Vq = —0.3. Monte-Carlo (1) uses (0), 
Monte-Carlo (2) uses 
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scattering length for V0=0.1 
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FIG. 5. The scattering length a for Vq = 0.1 as a function of the simulation time T at (3 
the exact value is a = —0.21. 
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FIG. 6. Probability distribution of path midpoints for Vq = 0.1 and o(r) as a function of the 
distance r. 
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Figure (1), the free action Figure (2), harmonic oscillator 




10 



10 100 10 100 

physical time (beta) physical time (beta) 

FIG. 7. The autocorrelation time tac for the free action. Figure (1): tac as a function of the 

physical time /3 for e = 0.5. Figure (2): the free action plus a weak harmonic oscillator potential 

(a; = 0.1). 
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l/e(T) for one and four nucleons 
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FIG. 8. The ratio of cross sections of the full problem and the effective potential l/e{T) as a 
function of the Langevin time T for nucleon and for four nucleons. 
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